#ifndef GEOMETRY_H
#define GEOMETRY_H

#include <math.h>

typedef struct { double x, y; } Point;
typedef struct { double a, b, c; } Line;
#if 0
typedef struct { Point O; double r; } Circle;
#else
typedef struct { Point O; double rsq; } Circle;
#endif

/* Returns the square of Euclidean distance between two points. */
double distance_squared(const Point P1, const Point P2)
{
    double dx = P1.x - P2.x;
    double dy = P1.y - P2.y;
    return dx * dx + dy * dy;
}

/* Returns the Euclidean distance between two points. */
double distance(const Point P1, const Point P2)
{
    double dx = P1.x - P2.x;
    double dy = P1.y - P2.y;
    return sqrt(dx * dx + dy * dy);
}

/* Returns the mid-point P of line segment P1P2 such that P1P / P1P2 = r. */
Point mid_point(const Point P1, const Point P2, double r)
{
    Point P;
    P.x = r*P2.x + (1-r)*P1.x;
    P.y = r*P2.y + (1-r)*P1.y;
    return P;
}

Line perpendicular_bisector(const Point P1, const Point P2)
{
    double x1 = P1.x, y1 = P1.y;
    double x2 = P2.x, y2 = P2.y;
    Line L;
    L.a = x2 - x1;
    L.b = y2 - y1;
    L.c = ((x1+x2)*(x1-x2)+(y1+y2)*(y1-y2))/2;
    return L;
}

Point circle_center(Circle C)
{
    return C.O;
}

double circle_radius(Circle C)
{
    return sqrt(C.rsq);
}

/* Returns a circle which is the locus of points P such that the ratio of 
 * distances from P to P1 and to P2 equals 'ratio'. 'ratio' must be positive
*  and not equal to 1.
 */
Circle apollonian_circle(const Point P1, const Point P2, double ratio)
{
    double x1 = P1.x, y1 = P1.y, x2 = P2.x, y2 = P2.y;
    double t = ratio * ratio;
    Circle C;
    C.rsq = ((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2))*t/((1-t)*(1-t));
    C.O.x = (x1-t*x2)/(1-t);
    C.O.y = (y1-t*y2)/(1-t);
    return C;
}

/* Intersects two lines. If there is a unique crossing point, stores the point
 * in P[0] and returns 1. Otherwise, returns 0.
 */
int intersect_lines(const Line L1, const Line L2, Point P[1])
{
    double a1 = L1.a, b1 = L1.b, c1 = L1.c;
    double a2 = L2.a, b2 = L2.b, c2 = L2.c;
    double det = a1*b2 - a2*b1;
    double p = b1*c2 - b2*c1;
    double q = a2*c1 - a1*c2;

    if (det == 0) /* parallel or overlapping */
        return 0;

    P[0].x = p / det;
    P[0].y = q / det;
    return 1;
}

/* Returns the perpendicular foot of a point to a line. */
Point perpendicular_foot(const Point P, const Line L)
{
    double a = L.a, b = L.b, c = L.c, x0 = P.x, y0 = P.y;
    double u = b*x0 - a*y0, v = a*a + b*b;
    Point M;
    M.x = (b*u-a*c)/v;
    M.y = (-a*u-b*c)/v;
    return M;
}

/* Intersects a line with a circle. If there are two (possibly overlapping)
 * crossing points, they are stored in P[0] and P[1] and the function returns
 * 2. Otherwise, returns 0.
 */
int intersect_line_circle(const Line L, const Circle C, Point P[2])
{
    double x0 = C.O.x, y0 = C.O.y;
    double a = L.a, b = L.b, c = L.c + a*x0 + b*y0;
    double delta = (a*a+b*b)*C.rsq - c*c;
    if (delta >= 0)
    {
        P[0].x = x0 + (-a*c+b*sqrt(delta))/(a*a+b*b);
        P[0].y = y0 + (-b*c-a*sqrt(delta))/(a*a+b*b);
        P[1].x = x0 + (-a*c-b*sqrt(delta))/(a*a+b*b);
        P[1].y = y0 + (-b*c+a*sqrt(delta))/(a*a+b*b);
        return 2;
    }
    return 0;
#if 0
    Point M = perpendicular_foot(circle_center(C), L);
    double d2 = distance_squared(circle_center(C), M);
    if (d2 <= C.rsq)
    {
        double u = sqrt((C.rsq - d2) / (L.a*L.a + L.b*L.b));
        P[0].x = M.x - L.b * u;
        P[0].y = M.y + L.a * u;
        P[1].x = M.x + L.b * u;
        P[1].y = M.y - L.a * u;
        return 2;
    }
    return 0;
#endif
}

/* Returns the radical line of two circles. */
Line radical_line(const Circle C1, const Circle C2)
{
    double x1 = C1.O.x, y1 = C1.O.y, r1sq = C1.rsq;
    double x2 = C2.O.x, y2 = C2.O.y, r2sq = C2.rsq;
    Line L;
    L.a = 2*(x1-x2);
    L.b = 2*(y1-y2);
    L.c = (r1sq-r2sq) - (x1+x2)*(x1-x2) - (y1+y2)*(y1-y2);
    return L;
}

/* Intersects two circles. If there are two (possibly overlapping) crossing
 * points, they are stored in P[0] and P[1] and the function returns 2. 
 * Otherwise, returns 0.
 */
int intersect_circles(const Circle C1, const Circle C2, Point P[2])
{
    return intersect_line_circle(radical_line(C1, C2), C1, P);
}

/* Returns the circumcenter of a triangle. If the triangle is degenerate,
 * the return value is undefined.
 */
Point circumcenter(const Point P[3])
{
    double x0 = P[0].x, y0 = P[0].y;
    double x1 = P[1].x, y1 = P[1].y;
    double x2 = P[2].x, y2 = P[2].y;
    double t = ((x1-x0)*(x2-x0)+(y1-y0)*(y2-y0)) / 
               ((x1-x0)*(y2-y0)-(x2-x0)*(y1-y0));
    Point O;
    O.x = ((x1+x2)+t*(y1-y2)) / 2;
    O.y = ((y1+y2)-t*(x1-x2)) / 2;
    return O;
}

#if 0
/* Returns a line that passes through the given point P and with the given
 * directional vector (u, v). */
Line make_ray(Point P, double u, double v)
{
    Line L;
    L.a = v;
    L.b = -u;
    L.c = u*P.y - v*P.x;
    return L;
}

/* Checks whether there exists a straight line that passes through three given
 * points.
 */
int colinear3(Point P1, Point P2, Point P3)
{
    double x1 = P1.x, y1 = P1.y;
    double x2 = P2.x, y2 = P2.y;
    double x3 = P3.x, y3 = P3.y;

}
#endif

#if 0
/* Returns a point with the given x and y coordinates. */
Point make_point(double x, double y)
{
    Point P;
    P.x = x;
    P.y = y;
    return P;
}

/* Makes a line that passes through two points. */
void make_line(Line L, const Point P1, const Point P2)
{
    L->a = P1->y - P2->y;
    L->b = P2->x - P1->x;
    L->c = -L->a * P1->x - L->b * P1->y;
}

/* Returns the square of Euclidean distance between two points. */
double distance_squared(const Point P1, const Point P2)
{
    double dx = P1->x - P2->x;
    double dy = P1->y - P2->y;
    return dx * dx + dy * dy;
}

/* Returns the sum of Euclidean distance between two sets of points. */
double sum_distance(const Point P1[], size_t n1, const Point P2[], size_t n2)
{
    double d = 0.0;
    size_t i, j;
    for (i = 0; i < n1; i++)
        for (j = 0; j < n2; j++)
            d += distance(P1[i], P2[j]);
    return d;
}

/* Finds the geometric median of four points, which is a point that minimizes
 * the total distance to the four points. Returns this minimum.
 */
double geometric_median(const Point P[4], Point M)
{
    double min_d = -1.0, d;
    int i;

    for (i = 0; i < 4; i++)
    {
        if ((d = sum_distance(&P[i], 1, P, 4)) < min_d || min_d < 0)
        {
            min_d = d;
            *M = *P[i];
        }
    }

    for (i = 0; i < 3; i++)
    {
        Line L1, L2;
        Point M2;
        make_line(L1, P[3], P[i]);
        make_line(L2, P[(i+1)%3], P[(i+2)%3]);
        if (intersect_lines(L1, L2, &M2) == 1)
            if ((d = sum_distance(&M2, 1, P, 4)) < min_d)
            {
                min_d = d;
                *M = *M2;
            }
    }
    
    return min_d;
}
#endif

#endif /* GEOMETRY_H */
